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Abstract. Generating accurate three dimensional planetary models and albedo maps is becoming 
increasingly more important as NASA plans more robotics missions to the Moon in the coming 
years. This paper describes a novel approach for separation of topography and albedo maps 
from orbital Lunar images. Our method uses an optimal Bayesian correlator to refine the stereo 
disparity map and generate a set of accurate digital elevation models (DEM). The albedo maps 
are obtained using a multi-image formation model that relies on the derived DEMs and the Lunar- 
Lambert reflectance model. The method is demonstrated on a set of high resolution scanned 
images from the Apollo era missions. 


1. Introduction 

High resolution, accurate topographic and albedo maps of planetary surfaces in general and Lunar 
surface in particular play an important role for the next NASA robotic missions. More specifically 
these maps are used in landing site selection, mission planing, planetary science discoveries and as 
educational resources. This paper describes a method for topographic and albedo maps reconstruc- 
tion from the Apollo era missions. The Apollo metric camera flown on an orbit at approximately 
100km above the Lunar surface was a calibrated wide field (75°) of view orbital mapping camera 
that photographed overlapping images (80%). The scans of these film images recently made avail- 
able [1, 2] capture the full dynamic range and resolution of the original film resulting in digital 
images of size 22,000 x 22,000 pixels representing a resolution of 10 m/pixel. Figure 1 shows the 
images of one Lunar orbit captured by the Apollo 15 mission. Our method for geometric stereo 



Figure 1 . Apollo Metric images from Orbit 33. 

reconstruction and photometric albedo reconstruction is illustrated in Figure 2. Each component of 
our system will be described in more detail in the following sections. 

2. Bundle Adjustment 

The Apollo-era satellite tracking network was highly inaccurate by today’s standards with errors 
estimated to be 2.04-km for satellite station positions and 0.002 degrees for pose estimates in a typical 
Apollo 15 image [3]. Such errors propagate through the stereo triangulation process, resulting in sys- 
tematic position errors and distortions in the resulting DEMs ( Figure 3). These errors are corrected 


*Carnegie Mellon University /NAS A Ames, ara.nefian@nasa.gov 

**NASA Ames, taemin.kim@nasa.gov, michael.broxton@nasa.gov, zach.moratto@nasa.gov. 



Orbital Image Pair 


Dense disparity maps ► DEM 



Figure 2. The overall system for albedo reconstruction. 



Figure 3. Bundle adjustment is illustrated here using a color- mapped, hill-shaded 
DEM mosaic from Apollo 15 Orbit 33 imagery, (a) Prior to bundle adjustment, 
large discontinuities exist between overlapping DEMs. (b) After bundle adjustment, 

DEM alignment errors are no longer visible. 

using bundle adjustment techniques. Our bundle adjustment solution uses SURF feature points [4]. 
Our bundle adjustment approach follows the method described in [5] and determines the best camera 
parameters that minimize the projection error given by e = ~ /(Cj, X&)) 2 where Ik are 

feature locations on the image plane, Cj are the camera parameters, and X & are the 3D positions 
associated with features /&. I{Cj,Xk) is an image formation model (i.e. forward projection) for a 
given camera and 3D point. The optimization of the cost function uses the Levenberg-Marquartd 
algorithm. Speed is improved by using sparse methods described in [6]. Outliers are rejected using 
the RANSAC method and trimmed to 1000 matches that are spread evenly across the images. To 
eliminate the gauge freedom inherent in this problem, we add two addition error metrics to this cost 
function to constrain the position and scale of the overall solution. First, e = — Cj) 2 

constrains camera parameters to stay close to their initial values. Second, a set of 3D ground control 
points are added to the error metric as e = J2 k (X k cp — X k ) 2 to constrain these points to known 
locations in the lunar coordinate frame. In the cost functions discussed above, errors are weighted 
by the inverse covariance of the measurement that gave rise to the constraint. Figure 3 shows a 
Lunar orbital DEM before and after the bundle adjustment processing. 

3. Dense Disparity Maps 

Apollo images are affected by two types of noise inherent to the scanning process: (1) the presence 
of film grain and (2) dust and lint particles. The former gives rise to noise in the DEM values that 
wash out real features, and the latter causes incorrect matches or hard to detect blemishes in the 
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DEM. Attenuating the effect of these scanning artifacts while simultaneously refining the integer 
disparity map to sub-pixel accuracy has become a critical goal of our system, and is necessary for 
processing real-world data sets such as the Apollo Metric Camera data. 

We investigated a large number of stereogrammmetric systems that can provide dense stereo 
matching from orbital imagery [7, 8, 9, 10, 11, 12]. A common technique in sub-pixel refinement is 
to fit a parabola to the correlation cost surface in the 8-connected neighborhood around the integer 
disparity estimate, and then use the parabola’s minimum as the sub-pixel disparity value. This 
method is easy to implement and fast to compute, but exhibits a problem known as pixel-locking: 
the sub-pixel disparities tend toward their integer estimates and can create noticeable ’’stair steps” 
on surfaces that should be smooth [12], [11]. One way of attenuating the pixel-locking effect is 
through the use of a symmetric cost function [8] for matching the “left” and “right” image blocks. 

To avoid the high computational complexity of these methods another class of approaches based on 
the Lucas-Kanade algorithm [13] proposes an asymmetric score where the disparity map is computed 
using the best matching score between the left image block and an optimally affine transformed 
block from the right image. For example, the sub-pixel refinement developed by Stein et. al. [12] 
lets /^(ra, n) and /^(i, j) be two corresponding pixels in the right and left image respectively, where 
i = m + d x , j = n + d y and d x , d y are the integer disparities. They develop a linear approximation 
based on the Taylor Series expansion around pixel (z, j) in the left image 
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where S x and S y are the local sub-pixel displacements. Let e(x,y) = In(x,y) — Il( i + S x ,j + S y ) 
and W be an image window centered around pixel (m, n). The local displacements are not constant 
across W and they vary according to: 
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The goal is to find the parameters ai, &i, ci, <22, 62, o 2 that minimize the cost function 
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where w(x,y) are a set of weights used to reject outliers. Note that the local displacements S x (i,j) 
and S y (i, j ) depend on the pixel positions within the window W. In fact, the values <21, 61, ci, <22, 62, c 2 
that minimize E can be seen as the parameters of an affine transformation that best transforms the 
right image window to match the reference (left) image window. 

The shortcoming of this method is directly related to the cost function that it is minimizing, 
which has a low tolerance to noise. Noise present in the image will easily dominate the result of the 
squared error function, giving rise to erroneous disparity information. Recently, several statistical 
approaches (e.g. [7]) have emerged to show how stochastic models can be used to attenuate the 
effects of noise. Our sub-pixel refinement technique [14] adopts some of these ideas, generalizing the 
earlier work by Stein et. al. [12] to a Bayesian framework that models both the data and image 
noise. 

In our approach the probability of a pixel in the right image is given by the following Bayesian 
model: 
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The first mixture component (z = 0) is a normal density function with mean Il(i + S x ,j + S y ) and 
variance Zf— : 
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The -p== factor in the variance of this component has the effect of a Gaussian smoothing window 
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over the patch. With this term in place, we are no longer looking for a single variance over the 
whole patch; instead we are assuming the variance increases with distance away from the center 
according to the inverted Gaussian, and are attempting to fit a global scale, a p . This provides 
formal justification for the standard Gaussian windowing kernel. 

The second mixture component (z = 1) in Equation 5 models the image noise using a normal 
density function with mean fi n and variance a n : 

(6) P(I R (m, n)\z = 1) = N(I R {m,n)\n n ,(j n ) 

Let I R (m,n) be a vector of all pixels values in a window W centered in pixel ( m,n ) in the right 
image. Then, 

(7) P(I R (m,n))= n P(Ir(x,v)) 
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The parameters A = {ai, fei , ci, C 2 , cr p , fi ni a n } that maximize the model likelihood in Equa- 

tion 7 are determined using the Expectation Maximization (EM) algorithm. Maximizing the model 
likelihood in Equation 7 is equivalent to maximizing the auxiliary function: 

Q ( 0 ) = 53 p ( fc |i*> A ‘)i°g p (i*>M A ) 
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Note that the M step calculations are similar to the equation used to determine the parameters 
ci, a 2 , 62 , C 2 in the method presented in [ 12 ], except here the fixed set of weights is replaced 
by the a posteriori probabilities computed in the E step. In this way, our approach can be seen as a 
generalization of the Lucas-Kanade method. The complete algorithm is summarized in the following 
steps: 

• Step 1 : Compute and the Ir(x , y) values using bilinear interpolation. 

Initialize the model parameters A. 

• Step 2: Compute iteratively the model parameters A using the EM algorithm (see [14] for 
details). 

• Step 3: Compute S x (i,j) and S y (i,j) using Equation 2 . 

• Step 4: Compute a new point (x f ,y f ) = (x,y) + (S x ,S y ) and the Ir(x' , y f ) values using 
bilinear interpolation. 

• Step 5: If the norm of (S x , S y ) vector falls below a fixed threshold the iterations converged. 
Otherwise, go to step 1 . 

Like the computation of the integer disparity maps, we adopt a multi-scale approach for sub-pixel 
refinement. At each level of the pyramid, the algorithm is initialized with the disparity determined 
in the previous lower resolution level of the pyramid. Figure 4 shows an example of a stereo image 
pair captured by the Apollo Metric Camera and used to generate a DEM of the Apollo 15 landing 
site. 


4. Photometric Reconstruction 

Each pixel of the Apollo Metric Camera images was formed by a combination of many factors, 
including albedo, terrain slope, exposure time, shadowing, and viewing and illumination angles. 
The goal of albedo reconstruction is to separate contributions of these factors. This is possible in 
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Figure 4. Apollo Metric Camera stereo pair showing Hadley Rille and the Apollo 
15 landing site: (left) left image, (middle) right image, (right) oblique view of the 
resulting DEM. 


part because of redundancy in the data; specifically, the same surface location is often observed in 
multiple overlapping images. 

To do the albedo reconstruction, we include all of the factors in a image formation model. Many 
of the parameters in this model such as digital terrain slopes, viewing angle, and sun ephemeris 
are known. To reconstruct albedo, we first model how the Metric Camera images were formed as a 
function of albedo, exposure time, illumination and viewing angles, and other factors. Then we can 
formulate the albedo inference problem as a least-squares solution that calculates the most likely 
albedo to produce the observed image data. 

Starting with the first images from the Apollo missions a large number of Lunar reflectance models 
were studied [15, 16, 17]. In this paper the reflectance is computed using the Lunar-Lambertian 
model [15, 18]. As shown in Figure 5, we define the following unit vectors: n is the local surface 
normal; 1 and v are directed at the locations of the Sun and the spacecraft, respectively, at the time 
when the image was captured. We further define the angles i separating n from 1 , e separating n 
from v, and the phase angle a separating 1 from v . The Lunar-Lambertian reflectance model is 
given by 


(9) 


F = AR = A 


(1 — L(a)) cos (i) + 2 L(a) 


cos (i) 

cos (i) + cos(e) 


where A is the intrinsic albedo and L(a) is a weighting factor between the Lunar and Lambertian 
reflectance models [19] that depends on the phase angle and surface properties. R is a photometric 
function that depends on the angles <a, i and e. The image formation model begins as follows. Let 



Figure 5. Illumination and viewing angles used by the Lunar-Lambertian re- 
flectance model. 
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Iij , Aij, Rij be the pixel value, albedo and R function at image location (i,j), and T be a variable 
proportional to the exposure time of the image. Then 

( 10 ) Rj = TAijRij. 

Note that the image formation model described in Equation 10 does not take into consideration the 
camera transfer function since the influence of the non-linearities of the camera transfer function 
plays a secondary role in the image formation model [19]. From Equation 10 it can be seen that 
when the observed pixel value, exposure time, and R value are known, the image formation model 
in Equation 10 provides a unique albedo value. However, these values are subject to errors arising 
from measurement (exposure time), scanning (image value) or stereo modeling errors (reflectance), 
resulting in imprecise albedo calculations. The method proposed here mitigates these errors by 
reconstructing the albedo of the Lunar surface from all the overlapping images, along with their 
corresponding exposure times and DEM information. The albedo reconstruction is formulated as 
the least squares problem that minimizes the following cost function Q: 

(11) Q = EE [(4 - AijT k Rfy 2 

k ij 

where super script k denotes the variables associated with the kth image and S^- is a shadow binary 
variable. Sfj = 1 when the pixel is in shadow and 0, otherwise. The weights w\- are chosen such that 
they have linearly decreasing values from the center of the image = 1) to the image boundaries 
(w!?j = 0). The choice of these weights insures that the reconstructed albedo mosaic is seamless. 
As shown by Equation 11 and illustrated in Figure 2 the steps of our photometric reconstruction 
method are the computation of the shadow and relief map followed by albedo reconstruction. These 
steps are described next. 



FIGURE 6 . Orbital image: (left) input image, (middle) binary shadow map with 
shadow regions shown in white, (right) DEM confidence map (brighter areas have 
higher estimated error). 


4.1. Shadow map computation. Discarding unreliable image pixels that are in shadow and for 
which the DEM and the reflectance models are unreliable plays an important role in accurate albedo 
estimation [20, 21]. Figure 6(left, middle) shows an input image together with its binary shadow 
map; shadowed areas are indicated in white. 

4.2. Relief map computation. The geodetically aligned local DEM determine multiple values 
for the same location on the Lunar surface. A simple average of the local DEM value determines 
the value used in computing the local slopes and the reflectance value. The average DEM has the 
following benefits for albedo reconstruction: 
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• It is essential to the computation of a coherent “ R map”, since each point of the Lunar 
surface must have a unique DEM value. 

• The statistical process produces more accurate terrain models by reducing the effect of 
random errors in local DEMs and without blurring the topographical features. Figure 7 
shows the R map of a subregion of the orbital image in Figure 6 before and after the DEM 
averaging and denoising process. It can be seen that the noise artifacts in the original DEM 
are reduced in the denoised DEM while the edges of the large crater and mountain regions 
are very well preserved. 

• The statistical parameters of the DEM values at each point are instrumental in building 
a confidence map of the Apollo coverage DEM. Figure 6 (right) shows the error confidence 
map for the orbital image illustrated in Figure 6(left). The values shown in this error map 
are the 0.05 x the variance values of the DEM expressed in meters. 

This step of the algorithm computes the values of the photometric function R described by 
Equation 9 corresponding to every pixel in the image. We denote the set of R values as the “R 
map” of the image. The accurate DEM calculation influences the R values through the effect of 
surface normals on the angles i and e that appear in Equation 9. 




Figure 7. R maps generated using (left) single local DEM and (right) denoised 
DEM derived from multiple overlapping local DEM. Our denoising approach pre- 
serves structure while reducing the artifacts shown in the insets. 



Figure 8. Albedo reconstruction: (left) R map, (middle) reconstructed albedo, 
(right) albedo confidence map (brighter areas have higher estimated error). 
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Figure 9. Albedo reconstruction of orbit 33 of the Apollo 15 mission. 


4.3. Albedo Reconstruction. The optimal albedo reconstruction [22] from multi view images and 
their corresponding DEM is formulated as a minimization problem of finding 

(12) {Aij,T k } = arg min Q 

Aij,T k 

for all pixels ij and images where Q is the cost function in Equation 11. An iterative solution to 
the above least square problem is given by the Gauss Newton updates described below. 

Step 1: Initialize the exposure time with the value provided in the image metadata. Ini- 
tialize the albedo map with the average value of the local albedo observed in all images. 

rk w k 

A . — V 
U uk rpk 

k U 

Step 2 : Re-estimate the albedo and refine the exposure time using 

i ,£ k (I? j -T>‘A ij R* j )T>'R* j S! j w$ i 

ij ij E k (T*Rk)*S*wk 

fk k 

Step 3: Compute the error cost function Q (Eqn. 11) for the re-estimated values of the 
albedo and exposure time. 

Convergence: If the convergence error between consecutive iterations falls below a fixed 
threshold then stop iterations and the re-estimated albedo is the optimal reconstructed 
albedo surface. Otherwise return to step 2. 

Figure 8 shows the R map, the albedo map and the albedo reconstruction error map, respectively, 
for the original orbital image in Figure 6. The albedo reconstruction error map is computed as 
the absolute difference between the original image 1^ and the reconstructed image T k AijR \y For 
display, the error values were multiplied by a factor 10 in Figure 8 (right). Figure 9 illustrates the 
reconstructed albedo for one orbit of the Apollo mission data overlayed over previous low resolution 
Clementine imagery. The Clementine mission images were captured under incidence and emission 
angles close to zero, therefore capturing images that describe the relative Lunar albedo. It can be 
seen that the reconstructed albedo de-emphasizes the brightness variations shown in the original 
imagery (Figure 1) between images and produces a seamless albedo mosaic. 

5. Conclusions 

This paper presents a novel approach for topographic and albedo maps generation from orbital 
imagery. The method for sub-pixel disparity maps uses a novel statistical formulation for optimally 
determining the stereo correspondence and reducing the effect of image noise. Our approach out- 
performs existing robust methods based on Lucas Kanade optical flow formulations at the cost of a 
higher computational complexity. The derived topographic maps are used to determine the albedo 
maps from an image formation model that incorporates the Lunar- Lambertian reflectance model. 


(13) 

(14) 

(15) 



The optimal values of albedo and exposure time are learned from multiple image views of the same 
area on Luna surface. Further research will be directed towards a joint estimation of the topgraphic 
and albedo information using shape from shading techniques specific for the Lunar reflectance model 
and scanned image properties. 
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